home *** CD-ROM | disk | FTP | other *** search
/ Perusromppu 2 / Perusromppu 2 - 40 suomenkielistä opetusohjelmaa.iso / perusr2 / lacom / lacom.for < prev    next >
Text File  |  1996-09-26  |  37KB  |  1,259 lines
  1.     PROGRAM LACOM
  2.  
  3. C Gas phase chemistry by Liisa Pirjola (box model)
  4. c for the Board of Education
  5. c solar intensity is an array of 79*57 called aflux3.dat 
  6. c zenith angle 30-86 astetta is enough in Finland
  7.  
  8.     IMPLICIT DOUBLE PRECISION ( A - Z )
  9.     INTEGER N,MAX,I,IW,IFAIL,I1,II,I7,I4,KK,I33,par
  10.     integer YEAR,day,mon,hour,pv,ku,tu,cl,nh,koodi,type1,nh1,lkm
  11.     PARAMETER (N=67, IW=(12+N)*N+50)
  12.     DOUBLE PRECISION TT,TEND,TOL
  13.     DOUBLE PRECISION W(IW),Y(N)
  14.     DOUBLE PRECISION J(79,57),ALFA
  15.     COMMON /INTENS/ J,ALFA
  16.     SAVE /INTENS/
  17.     DOUBLE PRECISION K(123),KF(16),H2O,O2,N2,H2,M
  18.     COMMON /LOPUT/ K,KF,H2O,O2,N2,H2,M
  19.     SAVE /LOPUT/
  20.     DOUBLE PRECISION V(67),EMI(67),VD(67),EMI1(67),wout
  21.     COMMON /NOP/ V,EMI,VD,EMI1,wout
  22.     SAVE /NOP/
  23.     DOUBLE PRECISION T
  24.     COMMON / TEMP /T
  25.     SAVE / TEMP /
  26.     EXTERNAL D02EAF,FCN
  27.  
  28. c simulation time (max), boundary layer height (h), latitude (ALAT), 
  29. c longitude (ALON), starting day (DAY), month (MON), hour (HOUR), 
  30. c albedo (ALB), uv radiation coefficient (G)
  31.     OPEN(17,FILE='INPUT1.DAT', STATUS = 'OLD')
  32.     READ(17,*) MAX,H,ALAT,ALON,YEAR,DAY,MON,HOUR,ALB,G,par
  33. c ticu0 is the starting time from the beginning of the year
  34.     ticu0=timecum(YEAR,day,mon,hour)
  35. c        write(6,*) ticu0
  36.     TEND=1.D-5
  37.     TENDAPU=TEND
  38.     TT = 0.
  39.     RAJA=1.
  40.  
  41. C initial concentrations (molecule cm-3), deposition velocities (m s-1)
  42. C  and emission rates (molecule cm-2 s-1)
  43.       OPEN(16,FILE = 'INPUT2.DAT', STATUS='OLD')
  44.       DO 15 I=1,54
  45.        READ(16,*) I33,Y(I),VD(I),EMI1(I)
  46. 15    CONTINUE
  47.       do 16 i=55,67
  48.     y(i)=0
  49.     vd(i)=0
  50.     emi1(i)=0
  51. 16    continue
  52.  
  53. c calls subroutine which calculates the zenith angle
  54.       CALL SOLAR(ALAT,ALON,YEAR,DAY,MON,HOUR)
  55.       open(11,file='zeni.dat',status ='old')
  56.  
  57. C reads intensity of the radiation into an array (79x57), the intensity
  58. c depends on the zenith angle and the wave length of the radiation
  59.       OPEN(12,FILE = 'AFLUX3.DAT',STATUS = 'OLD')
  60.       DO 300 KK=1,79
  61.       READ(12,*) (J(KK,II), II=1,57)
  62. 300   CONTINUE
  63.       CLOSE(12)
  64. c the intensity is multiplied by albedo and uv-coefficient 
  65.       if (alb.ge.0.46) then
  66.     do 301 kk=1,79
  67.     do 302 II=1,57
  68.     J(KK,II)=(1+alb)*G*J(KK,II)
  69. 302   continue
  70. 301   continue
  71.       else
  72.        do 305 kk=11,79
  73.        do 306 ii=1,57
  74.        j(kk,ii)=(1+alb)*G*J(kk,ii)
  75. 306   continue
  76. 305   continue
  77.       endif
  78.  
  79. c      write(6,*) '1+alb= ',1+alb,'uv coeff= ',g
  80. C open the restore files and save the initial values
  81.       OPEN(13,file='test1.res',status='unknown')
  82.       OPEN(14,file='test2.res',status='unknown')
  83.       WRITE(13,'(28e12.5)') TT,(Y(I1),I1=1,27)
  84.       write(14,'(28e12.5)') TT,(Y(I1),I1=28,54)
  85.  
  86.  
  87. C do 200 main-loop begins
  88.       DO 200 I=1,MAX
  89.       TICU=TICU0+TT/86400.
  90. c      write(6,*) 'ticu = ',ticu
  91.       OPEN(20,FILE='O3KRATE.DAT',STATUS='OLD')
  92.       OPEN(21,FILE='NO2.DAT',STATUS='OLD')
  93.       OPEN(22,FILE='NO3.DAT',STATUS='OLD')
  94.       OPEN(24,FILE='HNO3.DAT',STATUS='OLD')
  95.       OPEN(25,FILE='H2O2.DAT',STATUS='OLD')
  96.       OPEN(26,FILE='HCHO.DAT',STATUS='OLD')
  97.       OPEN(27,FILE='CH3CHO.DAT',STATUS='OLD')
  98.       OPEN(28,FILE='CHOCHO.DAT',STATUS='OLD')
  99.       OPEN(29,FILE='CH3OOH.DAT',STATUS='OLD')
  100.       OPEN(30,FILE='CH3COCHO.DAT',STATUS='OLD')
  101.        
  102.       read(11,*) pv,ku,tu,alfa
  103.       alfa=real(nint(alfa))
  104. c      write(6,*) 'alfa= ',alfa
  105.       OPEN(10,FILE='TEMP.DAT', STATUS='OLD')
  106.       OPEN(15,FILE='HUMI.DAT', STATUS='OLD')
  107.       
  108. c calculates the temperature in K
  109. 228   read(10,*) PV,KU,TU,TE2
  110. c      write(6,*) 'min tai max lämpötila ',te2
  111.       TICU2=TIMECUM(YEAR,PV,KU,TU)
  112.       IF (TICU2.LE.TICU) THEN
  113.     TICU1=TICU2
  114.     TE1=TE2
  115.     GOTO 228
  116.       ENDIF
  117.       close(10)
  118.       T=TE1+(TE2-TE1)/12.*(TICU-TICU1)
  119.       T=T+273.15
  120. c      WRITE(6,*) 'lämpötila Kelvineinä= ',T
  121.  
  122. c calculates H2O-concentration in molecule cm-3
  123. 229   read(15,*) PV,KU,TU,RH1,TE1
  124.       TICU1=TIMECUM(YEAR,PV,KU,TU)
  125.       IF (TICU1.LE.TICU) THEN
  126. c       write(6,*) 'ticu1 = ',ticu1
  127.        RH=0.01*RH1
  128.        TE=TE1+273.15
  129.        GOTO 229
  130.       ENDIF
  131. c      WRITE(6,*) 'RH= ',RH,'lämpötila = ',TE
  132.       psat = DEXP(77.344913-7235.4247/TE-8.2*DLOG(TE)+0.0057113*TE)
  133.       H2O = psat*RH/(TE*1.381d-17)
  134. c      WRITE(6,*) 'vesi = ',H2O
  135.       close(15)
  136.  
  137. c calculates the cloudiness coefficient with which the intensity is 
  138. c multiplied
  139.     open(18,file='clou.dat', status='old')
  140. 230   read(18,*) pv,ku,tu,type1,nh1
  141.     ticu1=timecum(YEAR,pv,ku,tu)
  142.     if (ticu1.le.ticu) then
  143.       cl=type1
  144.       nh=nh1
  145.       goto 230
  146.     endif
  147. c      write(6,*) 'tyyppi= ',cl,'määrä= ',nh
  148.       close(18)
  149.       CALL CLOUDINESS(CL,NH,ALFA,C1,koodi)
  150.       if (koodi.eq.1) goto 800 
  151.       do 303 kk=1,79
  152.       do 304 II=1,57
  153.       J(KK,II)=c1*J(KK,II)
  154. 304   continue
  155. 303   continue
  156. c      write(6,*) 'c1= ',c1
  157.  
  158. C calculates atmospheric oxygen, nitrogen and hydrogen values 
  159. c in molecule cm-3 
  160.       O2=0.209460*2.46e19*298./T
  161.       N2 =0.780840*2.46E19*298./T
  162.       H2 =0.5*2.46E13*298./T
  163.       M = O2+N2
  164.  
  165. c calls subroutines, which calculate deposition velocities (s-1), 
  166. c emission rates (molecule cm-2 s-1) and reaction rate constants
  167. c s-1 for   cm3 molecule-1 s-1 for   cm6 molecule -2 s-1 for
  168.       CALL DEPOSNOP(VD,H,V)
  169.       CALL EMISNOP(EMI1,H,EMI)
  170.       CALL RATECONS(T,H2O,O2,N2,H2,M,K)
  171.       if (par.eq.1) then 
  172.     CALL WASHOUT(T,H2O,WOUT)
  173.       else
  174.        wout=0
  175.       endif
  176. C calls subroutines, which give fotodissosiation rate coefficients (s-1)
  177.       CALL O3(ALFA,J,T,KF(1),KF(2))
  178.       CALL NO2(ALFA,J,KF(3))
  179.       CALL NO3(ALFA,J,KF(4),KF(5))
  180.       CALL N2O5(ALFA,J,T,KF(6))
  181.       CALL H2O2(ALFA,J,KF(7))
  182.       CALL HNO3(ALFA,J,KF(8))
  183.       CALL HCHO(ALFA,J,KF(9),KF(10))
  184.       CALL CH3CHO(ALFA,J,KF(11))
  185.       CALL CH3COC2H5(ALFA,KF(12))
  186.       CALL CH3COCOCH3(ALFA,KF(13))
  187.       CALL CHOCHO(ALFA,J,KF(14))
  188.       CALL CH3OOH(ALFA,J,KF(15))
  189.       CALL CH3COCHO(ALFA,J,KF(16))
  190.  
  191.       if (i.eq.1) then
  192.     lkm=66
  193.       else 
  194.     lkm=60
  195.       endif
  196.     DO 999 I7=1,lkm
  197.     TOL=1.e-4
  198.     IFAIL = 0
  199.     write(6,'(E10.4)') TT
  200. C calls D02EAF library routine
  201.     CALL D02EAF(TT,TEND,N,Y,TOL,FCN,W,IW,IFAIL)
  202.     DO 233 I4 = 1,54
  203. c        IF (Y(I4).LE.0.) Y(I4)=0.
  204.     IF (Y(I4).LE.1.d-5) Y(I4)=1.d-5
  205. 233     CONTINUE
  206.     TT = TEND
  207.     IF (TEND.LT.RAJA) THEN
  208.     TEND = 10.0*TEND
  209.     ELSE
  210.     TEND = TEND+60.
  211.     END IF
  212.       IF (MOD(NINT(TT-3601),3600).EQ.0.and.nint(tt).ne.1) THEN
  213.        WRITE(13,'(28e12.5)') TT-1,(Y(I1),I1=1,27)
  214.        write(14,'(28e12.5)') TT-1,(Y(I1),I1=28,54)
  215.       endif
  216. 999   continue
  217. 200   CONTINUE
  218.       CLOSE(10)
  219.       CLOSE(11)
  220.       CLOSE(15)
  221.       CLOSE(16)
  222. 800   STOP
  223.       END
  224.  
  225. C***********************************************************************
  226. C Routine for evaluating right hand sides of differential equations.
  227.       SUBROUTINE FCN(TT,Y,F)
  228.       IMPLICIT DOUBLE PRECISION ( A - Z)
  229. C N is the number of the differential equations
  230.       INTEGER N,i3
  231.       PARAMETER (N = 67)
  232. C TT is the time (s)
  233.       DOUBLE PRECISION TT,Y(N),F(N)
  234.       DOUBLE PRECISION K(123),KF(16),H2O,O2,N2,H2,M
  235.       COMMON /LOPUT/ K,KF,H2O,O2,N2,H2,M
  236.       DOUBLE PRECISION V(67),EMI(67),VD(67),EMI1(67),wout
  237.       COMMON /NOP/ V,EMI,VD,EMI1,wout
  238.       DOUBLE PRECISION APU(139)        
  239.       DOUBLE PRECISION T
  240.       COMMON / TEMP /T
  241.       
  242. c      do 256 i3=1,67
  243. c      if (y(i3).lt.1.e-28) y(i3)=1.e-28
  244. c256   continue
  245.  
  246. C kinetics
  247.       APU(1) = K(1)*Y(3)*O2
  248.       APU(2) = K(2)*Y(3)*Y(5)
  249.       APU(3) = K(3)*Y(2)*O2
  250.       APU(4) = K(4)*H2O*Y(2)
  251.       APU(5) = K(5)*Y(1)*Y(5)
  252.       APU(6) = K(6)*Y(1)*Y(6)
  253.       APU(7) = K(7)*Y(1)*Y(9)
  254.       APU(8) = K(8)*Y(1)*Y(10)
  255.       APU(9) = K(9)*Y(5)*Y(7)
  256.       APU(10) = K(10)*Y(5)*Y(10)
  257.       APU(11) = K(11)*Y(6)*Y(7)
  258.       APU(12) = K(12)*Y(6)*Y(7)
  259.       APU(13) = K(13)*Y(6)*Y(9)
  260.       APU(14) = K(14)*Y(7)*Y(12)
  261.       APU(15) = K(15)*2*Y(7)**2
  262.       APU(16) = K(16)*Y(8)*H2O
  263.       APU(17) = K(17)*Y(8)
  264.       APU(18) = K(18)*Y(9)*Y(10)
  265.       APU(19) = K(19)*Y(9)*Y(12)
  266.       APU(20) = K(20)*Y(9)*H2
  267.       APU(21) = K(21)*Y(9)*Y(13)
  268.       APU(22) = K(22)*2*Y(10)**2
  269.       APU(23) = K(23)*Y(10)*Y(7)
  270.       APU(24) = K(24)*Y(10)*Y(7)
  271.       APU(25) = K(25)*Y(9)*Y(14)
  272.       APU(26) = K(26)*Y(19)*Y(14)
  273.       APU(27) = K(27)*Y(9)*Y(17)
  274.       APU(28) = K(28)*Y(19)*Y(5)
  275.       APU(29) = K(29)*2*Y(19)**2
  276.       APU(30) = K(30)*2*Y(19)**2
  277.       APU(31) = K(31)*Y(10)*Y(19)
  278.       APU(32) = K(32)*Y(9)*Y(20)
  279.       APU(33) = K(33)*Y(7)*Y(20)
  280.       APU(34) = K(34)*Y(9)*Y(15)
  281.       APU(35) = K(35)*Y(21)*Y(9)
  282.       APU(36) = K(36)*Y(22)*Y(5)
  283.       APU(37) = K(37)*Y(23)
  284.       APU(38) = K(38)*Y(23)*O2
  285.       APU(39) = K(39)*Y(22)*Y(19)
  286.       APU(40) = K(40)*Y(9)*Y(24)
  287.       APU(41) = K(41)*Y(25)*Y(6)
  288.       APU(42) = K(42)*Y(26)
  289.       APU(43) = K(43)*Y(25)*Y(5)
  290.       APU(44) = K(44)*Y(19)*Y(25)
  291.       APU(45) = K(45)*2*Y(22)**2
  292.       APU(46) = K(46)*2*Y(25)**2
  293.       APU(47) = K(47)*Y(27)*Y(9)
  294.       APU(48) = K(48)*Y(28)*Y(5)
  295.       APU(49) = K(49)*Y(29)*O2
  296.       APU(50) = K(50)*Y(29)
  297.       APU(51) = K(51)*Y(28)*Y(19)
  298.       APU(52) = K(52)*Y(30)*Y(9)
  299.       APU(53) = K(53)*Y(31)*Y(5)
  300.       APU(54) = K(54)*Y(31)*Y(19)
  301.       APU(55) = K(55)*Y(34)*Y(9)
  302.       APU(56) = K(56)*Y(33)*Y(5)
  303.       APU(57) = K(57)*Y(33)*Y(19)
  304.       APU(58) = K(58)*Y(34)*Y(1)
  305.       APU(59) = K(59)*Y(35)*Y(1)
  306.       APU(60) = K(60)*Y(35)*Y(1)
  307.       APU(61) = K(61)*Y(35)*Y(9)
  308.       APU(62) = K(62)*Y(36)*Y(5)
  309.       APU(63) = K(63)*Y(36)*Y(19)
  310.       APU(64) = K(64)*Y(9)*Y(37)
  311.       APU(65) = K(65)*Y(47)*Y(5)
  312.       APU(66) = K(66)*Y(38)*Y(9)
  313.       APU(67) = K(67)*Y(39)*Y(5)
  314.       APU(68) = K(68)*Y(9)*Y(41)
  315.       APU(69) = K(69)*Y(9)*Y(40)
  316.       APU(70) = K(70)*Y(9)*Y(42)
  317.       APU(71) = K(71)*Y(43)*Y(5)
  318.       APU(72) = K(72)*Y(9)*Y(44)
  319.       APU(73) = K(73)*Y(45)*Y(5)
  320.       APU(74) = KF(1)*Y(1)
  321.       APU(75) = KF(2)*Y(1)
  322.       APU(76) = KF(3)*Y(6)
  323.       APU(77) = KF(4)*Y(7)
  324.       APU(78) = KF(5)*Y(7)
  325.       APU(79) = KF(6)*Y(8)
  326.       APU(80) = KF(7)*Y(12)
  327.       APU(81) = KF(8)*Y(13)
  328.       APU(82) = KF(9)*Y(20)
  329.       APU(83) = KF(10)*Y(20)
  330.       APU(84) = KF(11)*Y(24)
  331.       APU(85) = KF(12)*Y(30)
  332.       APU(86) = KF(13)*Y(32)
  333.       APU(87) = KF(14)*Y(41)
  334.       APU(88) = KF(15)*Y(46)
  335.       APU(89) = KF(16)*Y(40)
  336.       APU(90) = K(74)*Y(18)*O2
  337.       APU(91) = K(75)*Y(48)*O2
  338.       APU(92) = K(76)*Y(50)*O2
  339.       APU(93) = K(77)*Y(49)*H2O
  340.       APU(94) = K(78)*Y(42)*Y(1)
  341.       APU(95) = K(79)*Y(42)*Y(7)
  342.       APU(96) = K(80)*Y(51)*Y(9)
  343.       APU(97) = K(81)*Y(51)*Y(7)
  344.       APU(98) = K(82)*Y(51)*Y(1)
  345.       APU(99) = K(83)*Y(52)*Y(9)
  346.       APU(100) = K(84)*Y(52)*Y(7)
  347.       APU(101) = K(85)*Y(52)*Y(1)
  348.       APU(102) = K(86)*Y(53)*Y(9)
  349.       APU(103) = K(87)*Y(53)*Y(7)
  350.       APU(104) = K(88)*Y(53)*Y(1)
  351.       APU(105) = K(89)*Y(54)*Y(9)
  352.       APU(106) = K(90)*Y(54)*Y(7)
  353.       APU(107) = K(91)*Y(54)*Y(1)
  354.       APU(108) = K(92)*Y(55)*Y(9)
  355.       APU(109) = K(93)*Y(55)*Y(9)
  356.       APU(110) = K(94)*Y(57)
  357.       APU(111) = K(95)*Y(56)*Y(6)
  358.       APU(112) = K(96)*Y(56)*Y(1)
  359.       APU(113) = K(97)*Y(56)*O2
  360.       APU(114) = K(98)*Y(59)
  361.       APU(115) = K(99)*Y(61)*Y(6)
  362.       APU(116) = K(100)*Y(64)
  363.       APU(117) = K(101)*Y(59)
  364.       APU(118) = K(102)*Y(59)*Y(5)
  365.       APU(119) = K(103)*Y(58)*Y(1)
  366.       APU(120) = K(104)*Y(58)*Y(6)
  367.       APU(121) = K(105)*Y(58)*Y(6)
  368.       APU(122) = K(106)*Y(58)*O2
  369.       APU(123) = K(107)*Y(61)
  370.       APU(124) = K(108)*Y(61)*Y(5)
  371.       APU(125) = K(109)*Y(60)*Y(6)
  372.       APU(126) = K(110)*Y(60)*Y(1)
  373.       APU(127) = K(111)*Y(60)*O2
  374.       APU(128) = K(112)*Y(63)
  375.       APU(129) = K(113)*Y(63)*Y(6)
  376.       APU(130) = K(114)*Y(65)
  377.       APU(131) = K(115)*Y(60)
  378.       APU(132) = K(116)*Y(62)
  379.       APU(133) = K(117)*Y(62)
  380.       APU(134) = K(118)*Y(63)*Y(5)
  381.       APU(135) = K(119)*Y(59)*Y(1)
  382.       APU(136) = K(120)*Y(59)*Y(6)
  383.       APU(137) = K(121)*Y(66)
  384.       APU(138) = K(122)*Y(56)*O2
  385.       APU(139) = K(123)*Y(59)*O2
  386.  
  387. C Differential equations for concentrations
  388. c O3
  389.       F(1) = APU(1)-APU(5)-APU(6)-APU(7)-APU(8)-APU(58)-APU(59)
  390.      1       -APU(60)-APU(74)-APU(75)-APU(94)-APU(98)-APU(101)
  391.      2       -APU(104)-APU(107)-APU(112)-APU(119)-APU(126)
  392.      2       -APU(135)-V(1)*Y(1)
  393. c O(D)
  394.       F(2) = -APU(3)-APU(4)+APU(74)-V(2)*Y(2)
  395. c O(P)
  396.       F(3) = -APU(1)-APU(2)+APU(3)+APU(75)+APU(76)+APU(77)-V(3)*Y(3)
  397. c O2
  398.       F(4)=0
  399. c NO
  400.       F(5) = -APU(2)-APU(5)-APU(9)-APU(10)+APU(11)-APU(28)-APU(36)
  401.      1       -APU(43)-APU(48)-APU(53)-APU(56)-APU(62)-APU(65)-APU(67)
  402.      2       -APU(71)-APU(73)+APU(76)+APU(78)-apu(96)-apu(99)-apu(102)
  403.      3       -apu(105)+APU(111)-APU(118)+APU(120)+APU(121)-APU(124)
  404.      4       +APU(125)-APU(134)-V(5)*Y(5)+emi(5)
  405. c NO2
  406.       F(6) = APU(2)+APU(5)-APU(6)-APU(12)-APU(13)+2*APU(9)+APU(10)
  407.      1       +2*APU(15)+APU(17)+APU(24)+APU(28)+APU(36)-APU(41)+APU(42)
  408.      2       +APU(43)+APU(48)+APU(53)+APU(56)+APU(62)+APU(65)+APU(67)
  409.      3       +APU(71)+APU(73)-APU(76)+APU(77)+APU(79)+APU(81)+apu(96)
  410.      4       +apu(99)+apu(102)+apu(105)-APU(111)-APU(115)+APU(116)
  411.      5       +APU(118)-APU(120)-APU(121)+APU(124)-APU(125)-APU(129)
  412.      6       +APU(130)+APU(134)-APU(136)+APU(137)-V(6)*Y(6)+EMI(6)
  413. c NO3
  414.       F(7) = APU(6)-APU(9)-APU(11)-APU(12)-APU(14)-2*APU(15)+APU(17)
  415.      1       +APU(21)-APU(23)-APU(24)-APU(33)-APU(77)-APU(78)+APU(79)
  416.      2       -APU(95)-APU(97)-APU(100)-APU(103)-APU(106)-V(7)*Y(7)
  417. c N2O5
  418.       F(8) = APU(12)-APU(16)-APU(17)-APU(79)-V(8)*Y(8)
  419. c OH
  420.       F(9) = 2*APU(4)-APU(7)+APU(8)+APU(10)-APU(13)-APU(18)-APU(19)
  421.      1       -APU(20)-APU(21)+APU(24)-APU(25)-APU(27)-APU(32)-APU(34)
  422.      2       -APU(35)-APU(40)-APU(47)-APU(52)-APU(55)+0.19*APU(59)
  423.      3       -APU(61)-APU(64)-apu(66)-APU(68)-APU(69)-APU(70)-APU(72)
  424.      4       +2*APU(80)+APU(81)+APU(88)-APU(108)-APU(109)-APU(96)
  425.      5       -APU(102)-APU(105)-apu(99)-V(9)*Y(9)
  426. c HO2
  427.       F(10) = APU(7)-APU(8)-APU(10)+APU(14)-APU(18)+APU(19)-2*APU(22)
  428.      1       -APU(23)-APU(24)+APU(91)-APU(31)+apu(32)+APU(33)
  429.      2       +apu(34)+APU(38)+APU(49)+APU(51)+APU(53)+2*APU(54)
  430.      3       +APU(56)+0.12*APU(58)+0.29*APU(59)+0.12*APU(60)
  431.      4       +APU(62)+APU(63)+APU(65)+APU(67)+APU(71)+APU(73)
  432.      5       +2*APU(82)+APU(84)+APU(89)+APU(92)+apu(96)+apu(99)
  433.      6       +apu(102)+apu(105)-V(10)*Y(10)
  434. c H2SO4
  435.       F(11) = APU(93)-V(11)*Y(11)-wout*Y(11)
  436. c H2O2
  437.       F(12) = -APU(14)-APU(19)+APU(22)-APU(80)-V(12)*Y(12)-wout*Y(12)
  438. c HNO3
  439.       F(13) = APU(13)+APU(14)+2*APU(16)-APU(21)+APU(23)+APU(33)
  440.      1        -APU(81)-V(13)*Y(13)-wout*Y(13)
  441. c SO2
  442.       F(14) = -APU(25)-APU(26)+APU(121)+APU(131)+APU(138)+APU(139)
  443.      1         -V(14)*Y(14)+EMI(14)
  444. c       F(14) = 0
  445. c CO
  446.       F(15) = APU(32)+APU(33)-APU(34)+0.42*APU(58)+0.24*APU(59)
  447.      1        +0.42*APU(60)+2*APU(68)+APU(69)+APU(82)+APU(83)+APU(84)
  448.      2        +APU(87)+APU(89)-V(15)*Y(15)+emi(15)
  449. c CO2
  450.       F(16)=0
  451. c CH4
  452.       F(17)=0
  453. c CH3
  454.       F(18) = APU(27)-APU(90)+APU(37)+APU(43)+APU(44)+2*APU(46)
  455.      1        +APU(138)+APU(139)+APU(131)+APU(132)-V(18)*Y(18)
  456. c CH3O2
  457.       F(19) = -APU(26)+APU(90)-APU(28)-2*APU(29)-2*APU(30)-APU(31)
  458.      1        -APU(39)-APU(44)-APU(51)-APU(54)+0.43*APU(59)-APU(63)
  459.      2        +APU(84)-APU(57)-V(19)*Y(19)
  460. c HCHO
  461.       F(20) = APU(91)+APU(30)-APU(32)-APU(33)+APU(37)+APU(51)+APU(54)
  462.      1        +2*APU(56)+APU(58)+APU(59)+APU(62)+APU(63)
  463.      2        +APU(71)+APU(73)-APU(82)-APU(83)+APU(87)
  464.      3        -V(20)*Y(20)+EMI(20)
  465. c C2H6
  466.       F(21) = -APU(35)+EMI(21)-V(21)*Y(21)
  467. c C2H5O2
  468.       F(22) = APU(35)-APU(36)-APU(39)-2*APU(45)+APU(85)-V(22)*Y(22)
  469. c C2H5O
  470.       F(23) = APU(36)-APU(37)-APU(38)+APU(39)+2*APU(45)-V(23)*Y(23)
  471. c CH3CHO
  472.       F(24) = APU(38)-APU(40)+APU(50)+APU(60)+APU(62)-APU(84)+APU(63)
  473.      1        -V(24)*Y(24)+emi(24)
  474. c CH3COO2
  475.       F(25) = APU(40)-APU(41)+APU(42)-APU(43)-APU(44)-2*APU(46)+APU(69)
  476.      1        +APU(85)+2*APU(86)-V(25)*Y(25)
  477. c PAN
  478.       F(26) = APU(41)-APU(42)-V(26)*Y(26)
  479. c nC4H10
  480.       F(27) = -APU(47)+EMI(27)-V(27)*Y(27)
  481. c C4H9O2
  482.       F(28) = APU(47)-APU(48)-APU(51)-V(28)*Y(28)
  483. c C4H9O
  484.       F(29) = APU(48)-APU(49)-APU(50)+APU(51)-V(29)*Y(29)
  485. c CH3COC2H5
  486.       F(30) = APU(49)-APU(52)-APU(85)-V(30)*Y(30)+emi(30)
  487. c CH3COCHO2CH3
  488.       F(31) = APU(52)-APU(53)-APU(54)-V(31)*Y(31)
  489. c CH3COCOCH3
  490.       F(32) = APU(53)+APU(54)-APU(86)-V(32)*Y(32)
  491. c CH2O2CH2OH
  492.       F(33) = APU(55)-APU(56)-APU(57)-V(33)*Y(33)
  493. c C2H4
  494.       F(34) = -APU(55)-APU(58)+EMI(34)-v(34)*Y(34)
  495. c C3H6
  496.       F(35) = -APU(59)-APU(60)-APU(61)+EMI(35)-V(35)*Y(35)
  497. c CH3CHO2CH2OH
  498.       F(36) = APU(61)-APU(62)-APU(63)-V(36)*Y(36)
  499. c o-ksyleeni
  500.       F(37) = -APU(64)+EMI(37)-V(37)*Y(37)
  501. c CHOCH=CHCOCH3
  502.       F(38) = APU(65)-APU(66)-V(38)*Y(38)
  503. c CH3COCHOH-CHO2CHO
  504.       F(39) = APU(66)-APU(67)-V(39)*Y(39)
  505. c CH3COCHO
  506.       F(40) = APU(73)+APU(65)-APU(69)-APU(89)+APU(67)-V(40)*Y(40)
  507. c CHOCHO
  508.       F(41) = APU(67)-APU(68)-APU(87)-V(41)*Y(41)
  509. c C5H8
  510.       F(42) = -APU(70)-APU(94)-APU(95)+EMI(42)-V(42)*Y(42)
  511. c OHC5H8O2
  512.       F(43) = APU(70)-APU(71)-V(43)*Y(43)
  513. c CH3COCH=CH2
  514.       F(44) = APU(71)-APU(72)-V(44)*Y(44)
  515. c OHCH3COCHCH2O2
  516.       F(45) = APU(72)-APU(73)-V(45)*Y(45)
  517. c CH3OOH
  518.       F(46) = APU(31)-APU(88)-V(46)*Y(46)-wout*Y(46)
  519. c
  520.       F(47) = APU(64)-APU(65)
  521. c CH3O
  522.       F(48) = APU(26)+APU(28)+0.05*APU(59)+2*APU(29)+APU(39)+APU(44)
  523.      1        +APU(57)+APU(63)+apu(88)-APU(91)-V(48)*Y(48)
  524. C SO3
  525.       F(49) = APU(26)+APU(92)-APU(93)+APU(132)-V(49)*Y(49)
  526. C HSO3
  527.       F(50) = APU(25)-APU(92)-V(50)*Y(50)
  528. C ALFAPINENE
  529.       F(51) = -APU(96)-APU(97)-APU(98)+EMI(51)-V(51)*Y(51)
  530. C BEETAPINENE
  531.       F(52) = -APU(99)-APU(100)-APU(101)+EMI(52)-V(52)*Y(52)
  532. C 3-KARENE
  533.       F(53) = -APU(102)-APU(103)-APU(104)+EMI(53)-V(53)*Y(53)
  534. C D-LIMONENE
  535.       F(54) = -APU(105)-APU(106)-APU(107)+EMI(54)-V(54)*Y(54)
  536. C DMS
  537. c      F(55) = -APU(108)-APU(109)+EMI(55)-V(55)*Y(55)
  538.        f(55)=0
  539. C CH3S
  540. c      F(56) = APU(108)-APU(111)-APU(112)-APU(113)+APU(114)-APU(138)
  541. c     1       -V(56)*Y(56)
  542.       f(56)=0
  543. C CH3S(OH)CH3
  544. c      F(57) = APU(109)-APU(110)-V(57)*Y(57)
  545.       f(57)=0
  546. C CH3SO
  547. c      F(58) = APU(110)+APU(111)+APU(112)+APU(118)-APU(119)-APU(120)
  548. c     1        -APU(121)-APU(122)+APU(123)+APU(135)-V(58)*Y(58)
  549.       f(58)=0
  550. C CH3SOO
  551. c      F(59) = APU(113)-APU(114)-APU(117)-APU(118)-APU(135)-APU(136)
  552. c     1        +APU(137)-APU(139)-V(59)*Y(59)
  553.       f(59)=0
  554. C CH3SO2
  555. c      F(60) = APU(117)+APU(119)+APU(120)+APU(124)-APU(125)-APU(126)
  556. c     1        -APU(127)+APU(128)-V(60)*Y(60)
  557.       f(60)=0
  558. C CH3S(O)O2
  559. c     F(61) = -APU(115)+APU(116)+APU(122)-APU(123)-APU(124)-V(61)*Y(61)
  560.       f(61)=0
  561. C CH3SO3
  562. c      F(62) = APU(125)+APU(126)-APU(132)-APU(133)+APU(134)-V(62)*Y(62)
  563.       f(62)=0
  564. C CH3S(O)2O2
  565. c      F(63) = APU(127)-APU(128)-APU(129)+APU(130)-APU(134)-V(63)*Y(63)
  566.       f(63)=0
  567. C CH3S(O)O2NO2
  568. c      F(64) = APU(115)-APU(116)-V(64)*Y(64)
  569.       f(64)=0
  570. C CH2S(O)2O2NO2
  571. c      F(65) = APU(129)-APU(130)-V(65)*Y(65)
  572.       f(65)=0
  573. C CH3SOONO2
  574. c      F(66) = APU(136)-APU(137)-V(66)*Y(66)
  575.       f(66)=0
  576. C CH3SO3H (MSA)
  577. c      F(67) = APU(133)-V(67)*Y(67)
  578.       f(67)=0
  579. c      do 257 i3=1,54      
  580. c      if (f(i3).gt.1.e20) write(6,*) f(i3)
  581. c257   continue
  582.       RETURN
  583.       END
  584. C
  585.  
  586. C subroutines for the rate constants of the photodissociation reactions
  587.       SUBROUTINE O3(ALFA,J,T,K1,K2)
  588.       IMPLICIT DOUBLE PRECISION (A - Z)
  589.       INTEGER I,L
  590.       DOUBLE PRECISION J(79,57),ALFA
  591.       DOUBLE PRECISION T
  592.       K1=0.
  593.       DL=5.
  594.       VA=1.E-20*1.E14
  595.       T1=T-230.
  596.       A=0.332+2.565E-4*T1+1.152E-5*T1**2+2.313E-8*T1**3
  597.       B=-0.575+5.59E-3*T1-1.439E-5*T1**2+3.27E-8*T1**3
  598.       C=0.466+8.883E-4*T1-3.546E-5*T1**2-3.519E-7*T1**3
  599.       LO=308.20+4.4871E-2*T1+6.9380E-5*T1*2-2.5452E-6*T1**3
  600.       DO 500 I=1,7
  601.       READ(20,*) AL,CR
  602.       L=NINT(AL)
  603.       IF (L.LE.300) THEN
  604.     FI=0.9
  605.       ELSE
  606.     FI=A*ATAN(B*(REAL(L)-LO))+C
  607.       END IF
  608.       IF (FI.GT.0.9) FI=0.9
  609.       IF (FI.LT.0.) FI=0.
  610.       IF (ALFA.GT.86.) THEN
  611.        K1=K1+CR*FI*J((L-290)/5+1,57)*DL*VA
  612.        K1 = -K1/4*(ALFA-90.)
  613.       else
  614.        K1=K1+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA 
  615.       endif
  616. 500   CONTINUE
  617.       CLOSE(20)
  618.       IF (ALFA.GT.89.9) THEN
  619.        KULMA=89.9
  620.        K2=1.23E-3*DEXP(-0.6/COS(KULMA*3.14159265/180.))
  621.       ELSE 
  622.        K2=1.23E-3*DEXP(-0.6/COS(ALFA*3.14159265/180.))
  623.       ENDIF
  624.       if (alfa.ge.90.) k2=0.
  625.       RETURN
  626.       END
  627. C
  628.       SUBROUTINE NO2(ALFA,J,K)
  629.       IMPLICIT DOUBLE PRECISION (A - Z)
  630.       INTEGER I,L
  631.       DOUBLE PRECISION J(79,57),ALFA
  632.       K=0.
  633.       DL=5.
  634.       VA=1.E-20*1.E14
  635.       DO 500 I=1,25
  636.       READ(21,*) L,CR,FI
  637.       IF (ALFA.GT.86.) THEN
  638.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  639.        K = -K/4*(ALFA-90.)
  640.       ELSE
  641.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  642.       ENDIF
  643.   500 CONTINUE
  644.       CLOSE(21)
  645.       RETURN
  646.       END
  647. C
  648.       SUBROUTINE NO3(ALFA,J,K1,K2)
  649.       IMPLICIT DOUBLE PRECISION (A - Z)
  650.       INTEGER I,L
  651.       DOUBLE PRECISION J(79,57),K1,K2
  652.       K1=0.
  653.       K2=0.
  654.       DL=5.
  655.       VA=1.E-19*1.E14
  656.       DO 500 I=1,32
  657.       READ(22,*) AL,TULO1,TULO2
  658.       L=NINT(AL)
  659.       IF (ALFA.GT.86.) THEN
  660.        K1=K1+TULO1*J((L-290)/5+1,57)*DL*VA
  661.        K1 = -K1/4*(ALFA-90.)
  662.        K2=K2+TULO2*J((L-290)/5+1,57)*DL*VA
  663.        K2 = -K2/4*(ALFA-90.)
  664.       ELSE
  665.        K1=K1+TULO1*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  666.        K2=K2+TULO2*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  667.       ENDIF
  668.   500 CONTINUE
  669.       CLOSE(22)
  670.       RETURN
  671.       END
  672. C
  673.       SUBROUTINE N2O5(ALFA,J,T,K)
  674.       IMPLICIT DOUBLE PRECISION (A - Z)
  675.       INTEGER I,L
  676.       DOUBLE PRECISION J(79,57),ALFA
  677.       K=0.
  678.       DL=5.
  679.       L=290
  680.       FI=1.
  681.       VA=1.E-20*1.E14
  682.       DO 500 I=1,13
  683.       CR=DEXP(2.735+(4728-17.13*L)/T)
  684.       IF (ALFA.GT.86.) THEN
  685.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  686.        K = -K/4*(ALFA-90.)
  687.       ELSE
  688.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  689.       ENDIF
  690.       L=L+5
  691.   500 CONTINUE
  692.       RETURN
  693.       END
  694. C
  695.       SUBROUTINE HNO3(ALFA,J,K)
  696.       IMPLICIT DOUBLE PRECISION (A - Z)
  697.       INTEGER I,L
  698.       DOUBLE PRECISION J(79,57),ALFA
  699.       K=0.
  700.       DL=5.
  701.       VA=1.E-20*1.E14
  702.       DO 500 I=1,7
  703.       READ(24,*) L,CR,FI
  704.       IF (ALFA.GT.86.) THEN
  705.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  706.        K = -K/4*(ALFA-90.)
  707.       ELSE
  708.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  709.       ENDIF
  710.   500 CONTINUE
  711.       CLOSE(24)
  712.       RETURN
  713.       END
  714. C
  715.       SUBROUTINE H2O2(ALFA,J,K)
  716.       IMPLICIT DOUBLE PRECISION (A - Z)
  717.       INTEGER I,L
  718.       DOUBLE PRECISION J(79,57),ALFA
  719.       K=0.
  720.       DL=5.
  721.       VA=1.E-20*1.E14
  722.       DO 500 I=1,13
  723.       READ(25,*) L,CR,FI
  724.       IF (ALFA.GT.86.) THEN
  725.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  726.        K = -K/4*(ALFA-90.)
  727.       ELSE
  728.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  729.       ENDIF
  730.   500 CONTINUE
  731.       CLOSE(25)
  732.       RETURN
  733.       END
  734. C
  735.       SUBROUTINE HCHO(ALFA,J,K1,K2)
  736.       IMPLICIT DOUBLE PRECISION (A - Z)
  737.       INTEGER I,L
  738.       DOUBLE PRECISION J(79,57),ALFA
  739.       K1=0.
  740.       K2=0.
  741.       DL=5.
  742.       VA=1.E-20*1.E14
  743.       DO 500 I=1,15
  744.       READ(26,*) AL,CR,FI1,FI2
  745.       L=NINT(AL)
  746.       IF (ALFA.GT.86.) THEN
  747.        K1=K1+CR*FI1*J((L-290)/5+1,57)*DL*VA
  748.        K2=K2+CR*FI2*J((L-290)/5+1,57)*DL*VA
  749.        K1 = -K1/4*(ALFA-90.)
  750.        K2 = -K2/4*(ALFA-90.)
  751.       ELSE
  752.        K1=K1+CR*FI1*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  753.        K2=K2+CR*FI2*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  754.       ENDIF
  755.   500 CONTINUE
  756.       CLOSE(26)
  757.       RETURN
  758.       END
  759. C
  760.       SUBROUTINE CH3CHO(ALFA,J,K)
  761.       IMPLICIT DOUBLE PRECISION (A - Z)
  762.       INTEGER I,L
  763.       DOUBLE PRECISION J(79,57),ALFA
  764.       K=0.
  765.       DL=5.
  766.       VA=1.E-20*1.E14
  767.       DO 500 I=1,8
  768.       READ(27,*) L,CR,FI
  769.       IF (ALFA.GT.86.) THEN
  770.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  771.        K = -K/4*(ALFA-90.)
  772.       ELSE
  773.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  774.       ENDIF
  775.   500 CONTINUE
  776.       CLOSE(27)
  777.       RETURN
  778.       END
  779. C
  780.       SUBROUTINE CHOCHO(ALFA,J,K)
  781.       IMPLICIT DOUBLE PRECISION (A - Z)
  782.       INTEGER I,L
  783.       DOUBLE PRECISION J(79,57),ALFA
  784.       K=0.
  785.       DL=5.
  786.       VA=1.E-20*1.E14
  787.       DO 500 I=1,35
  788.       READ(28,*) AL,CR
  789.       L=NINT(AL)
  790.       IF (L.LT.325) THEN
  791.     FI=0.4
  792.       ELSE
  793.     FI=0.029
  794.       ENDIF
  795.       IF (ALFA.GT.86.) THEN
  796.     K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  797.     K = -K/4*(ALFA-90.)
  798.       ELSE
  799.     K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  800.       ENDIF
  801. 500   CONTINUE
  802.       CLOSE(28)
  803.       RETURN
  804.       END
  805. C
  806.       SUBROUTINE CH3OOH(ALFA,J,K)
  807.       IMPLICIT DOUBLE PRECISION (A - Z)
  808.       INTEGER I,L
  809.       DOUBLE PRECISION J(79,57),ALFA
  810.       K=0.
  811.       DL=5.
  812.       VA=1.E-20*1.E14
  813.       DO 500 I=1,13
  814.       READ(29,*) AL,CR,FI
  815.       L=NINT(AL)
  816.       IF (ALFA.GT.86.) THEN
  817.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  818.        K = -K/4*(ALFA-90.)
  819.       ELSE
  820.     K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  821.       ENDIF
  822.   500 CONTINUE
  823.       CLOSE(29)
  824.       RETURN
  825.       END
  826. C
  827.       SUBROUTINE CH3COCHO(ALFA,J,K)
  828.       IMPLICIT DOUBLE PRECISION (A - Z)
  829.       INTEGER I,L
  830.       DOUBLE PRECISION J(79,57),ALFA
  831.       K=0.
  832.       DL=5.
  833.       VA=1.E-20*1.E14
  834.       FI=0.107
  835.       DO 500 I=1,29
  836.       READ(30,*) AL,CR
  837.       L=NINT(AL)
  838.       IF (ALFA.GT.86.) THEN
  839.        K=K+CR*FI*J((L-290)/5+1,57)*DL*VA
  840.        K = -K/4*(ALFA-90.)
  841.       ELSE
  842.        K=K+CR*FI*J((L-290)/5+1,NINT(ALFA)-29)*DL*VA
  843.       ENDIF
  844.   500 CONTINUE
  845.       CLOSE(30)
  846.       RETURN
  847.       END
  848. C
  849.       SUBROUTINE CH3COC2H5(ALFA,K)
  850.       IMPLICIT DOUBLE PRECISION (A - Z)
  851.       DOUBLE PRECISION ALFA
  852.       A=2.43E-5
  853.       B=0.877
  854.       IF (ALFA.GT.89.9) THEN
  855.        K=A*DEXP(-B/COS(89.9*3.14159265/180.))
  856.       ELSE
  857.        K=A*DEXP(-B/COS(ALFA*3.14159265/180.))
  858.       ENDIF
  859.       if (alfa.ge.90.) k=0.
  860.       RETURN
  861.       END
  862. C
  863.       SUBROUTINE CH3COCOCH3(ALFA,K)
  864.       IMPLICIT DOUBLE PRECISION (A - Z)
  865.       DOUBLE PRECISION ALFA
  866.       A=2.7E-4
  867.       B=0.79
  868.       IF (ALFA.GT.89.9) THEN
  869.        K=A*DEXP(-B/COS(89.9*3.14159265/180.))
  870.       ELSE
  871.        K=A*DEXP(-B/COS(ALFA*3.14159265/180.))
  872.       ENDIF
  873.       if (alfa.ge.90.) k=0.
  874.       RETURN
  875.       END
  876.  
  877. C Rate constants
  878.       SUBROUTINE RATECONS(T,H2O,O2,N2,H2,M,K)
  879.       IMPLICIT DOUBLE PRECISION (A - Z)
  880.       DOUBLE PRECISION T,K(123),KF(16),O2,N2,M
  881.  
  882.     KoO2=6.2E-34*(T/300.)**(-2.8)*O2
  883.     KoN2=5.6E-34*(T/300.)**(-2.8)*N2
  884.     KINF=2.8E-12
  885.     Fc=DEXP(-T/696.)
  886.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  887.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  888.       K(1) = KO2+KN2
  889.     KoO2=8.6E-32*(T/300.)**(-1.8)*O2
  890.     KoN2=1.0E-31*(T/300.)**(-1.6)*N2
  891.     KINF=3.0E-11*(T/300.)**0.3
  892.     Fc=DEXP(-T/1850.)
  893.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  894.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  895.       K(2) = KO2+KN2
  896.       K(3) = 3.2E-11*DEXP(67./T)
  897.       K(4) = 2.2E-10
  898.       K(5) = 1.8E-12*DEXP(-1370./T)
  899.       K(6) = 1.2E-13*DEXP(-2450./T)
  900.       K(7) = 1.9E-12*DEXP(-1000./T)
  901.       K(8) = 1.4E-14*DEXP(-600./T)
  902.       K(9) = 1.8E-11*DEXP(110./T)
  903.       K(10) = 3.7E-12*DEXP(240./T)
  904.       K(11) = 2.3E-12*DEXP(-1000./T)
  905.     KoN2=2.7E-30*(T/300.)**(-3.4)*N2
  906.     KINF=2.0E-12*(T/300.)**0.2
  907.     Fc = DEXP(-T/250.)+DEXP(-1050./T)
  908.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  909.       K(12) = KN2
  910.     KoO2=2.2E-30*(T/300.)**(-2.9)*O2
  911.     KoN2=2.6E-30*(T/300.)**(-2.9)*N2
  912.     KINF=5.2E-11
  913.     Fc=DEXP(-T/353.)
  914.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  915.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  916.       K(13) = KO2+KN2
  917.       K(13) = KN2
  918.       K(14) = 4.1E-16
  919.       K(15) = 8.5E-13*DEXP(-2450./T)
  920.       K(16) = 1.3E-21
  921.     KoN2=2.2E-3*(T/300.)**(-4.4)*DEXP(-11080./T)*N2
  922.     KINF=9.7E14*(T/300.)**0.1*DEXP(-11080./T)
  923.     Fc=DEXP(-T/250.)+DEXP(-1050./T)
  924.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  925.       K(17) = KN2
  926.       K(18) = 4.8E-11*DEXP(250./T)
  927.       K(19) = 2.9E-12*DEXP(-160./T)
  928.       K(20) = 7.7E-12*DEXP(-2100./T)
  929.       K(21) = 9.4E-15*DEXP(778./T)
  930.       K(22) = 1.6E-12
  931.       K(23) = 4.3E-12
  932.       K(24) = 4.3E-12
  933.     KoO2=5.0E-31*(T/300.)**(-3.3)*O2
  934.     KoN2=5.0E-31*(T/300.)**(-3.3)*N2
  935.     KINF=2.E-12
  936.     Fc=DEXP(-T/380)
  937.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  938.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  939.       K(25) = KO2+KN2
  940.       K(25) = KO2
  941.       K(26) = 4.0E-17
  942.       K(27) = 3.9E-12*DEXP(-1885./T)
  943.       K(28) = 4.2E-12*DEXP(180./T)
  944.       K(29) = 1.1E-13*DEXP(365./T)
  945.       K(30) = 1.1E-13*DEXP(365./T)
  946.       K(31) = 3.8E-13*DEXP(780./T)
  947.       K(32) = 8.8E-12*DEXP(25./T)
  948.       K(33) = 5.8E-16
  949.       K(34) = 2.4E-13
  950.       K(35) = 7.8E-12*DEXP(-1020./T)
  951.       K(36) = 8.9e-12
  952.       K(37) = 33.0
  953.       K(38) = 9.5E-15
  954.       K(39) = 2.5E-14
  955.       K(40) = 5.6E-12*DEXP(310./T)
  956.     KoN2=2.7E-28*(T/300.)**(-7.1)*N2
  957.     KINF=1.2E-11*(T/300.)**(-0.9)
  958.     Fc=0.3
  959.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  960.       K(41) = KN2
  961.     KoN2=4.9E-3*DEXP(-12100./T)*N2
  962.     KoN2=4.9E-3*DEXP(-12100./T)*M
  963.     KINF=4.0E16*DEXP(-13600./T)
  964.     Fc=0.3
  965.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  966.       K(42) = KN2
  967.       K(43) = 2.0E-11
  968.       K(44) = 5.5E-12
  969.       K(45) = 9.8E-14*DEXP(-110./T)
  970.       K(46) = 2.8E-12*DEXP(530./T)
  971.       K(47) = 1.4E-11*DEXP(-559./T)
  972.       K(48) = 3.0E-12
  973.       K(49) = 2.1E-16
  974.       K(50) = 1.2E3
  975.       K(51) = 2.5E-14
  976.       K(52) = 8.8E-13
  977.       K(53) = 3.1E-13
  978.       K(54) = 2.5E-14
  979.     KoO2=9.5E-29*(T/300.)**(-3.1)*O2
  980.     KoN2=9.5E-29*(T/300.)**(-3.1)*N2
  981.     KINF=9.E-12
  982.     Fc=DEXP(-T/840.)
  983.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  984.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  985.       K(55) = KO2+KN2
  986.       K(55) = KN2
  987.       K(56) = 3.1E-13
  988.       K(57) = 2.5E-14
  989.       K(58) = 1.2E-14*DEXP(-2630./T)
  990.       K(59) = 6.5E-15*DEXP(-1880./T)
  991.       K(60) = 6.5E-15*DEXP(-1880./T)
  992.     KoO2=8.E-27*(T/300.)**(-3.5)*O2
  993.     KoN2=8.E-27*(T/300.)**(-3.5)*N2
  994.     KINF=3.0E-11
  995.     Fc=DEXP(-T/433.)
  996.     KO2=(KoO2/(1.+KoO2/KINF))*Fc**(1/(1+(DLOG10(KoO2/KINF))**2))
  997.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  998.       K(61) = KO2+KN2
  999.       K(62) = 3.1E-13
  1000.       K(63) = 2.5E-14
  1001.       K(64) = 1.1E-11
  1002.       K(65) = 3.1E-13
  1003.       K(66) = 2.0E-11
  1004.       K(67) = 3.1E-13
  1005.       K(68) = 1.1E-11
  1006.       K(69) = 1.7E-11
  1007.       K(70) = 2.54E-11*DEXP(410./T)
  1008.       K(71) = 3.0E-13
  1009.       K(72) = 2.0E-11
  1010.       K(73) = 3.0E-13
  1011.     KoN2=1.0E-30*(T/300.)**(-3.3)*N2
  1012.     KINF=2.2E-12*(T/300.)**1.0
  1013.     Fc=0.27
  1014.     KN2=(KoN2/(1.+KoN2/KINF))*Fc**(1/(1+(DLOG10(KoN2/KINF))**2))
  1015.       K(74) = KN2
  1016.       K(75) = 1.9E-15
  1017.       K(76) = 4.4E-13
  1018.       K(77) = 9.1E-13
  1019.       K(78) = 1.2e-17
  1020.       K(79) = 3.2e-13
  1021.       K(80) = 0.98e-11*dexp(500./T)
  1022.       K(81) = 5.8e-12
  1023.       K(82) = 4.6e-15*dexp(-1170./T)
  1024.       K(83) = 7.95e-11
  1025.       K(84) = 2.36e-12
  1026.       K(85) = 2.1e-17
  1027.       K(86) = 1.6e-11*dexp(500./T)
  1028.       K(87) = 10.1e-12
  1029.       K(88) = 6.5e-15*dexp(-1170./T)
  1030.       K(89) = 16.9e-11
  1031.       K(90) = 1.4e-11
  1032.       K(91) = 6.4e-16
  1033.       K(92) = 1.1D-11*dexp(-240./T)
  1034.       K(93) = 1.2D-12
  1035.       K(94) = 0.5
  1036.       K(95) = 5.6D-11
  1037.       K(96) = 5.4D-12
  1038.       K(97) = 6.1D-19
  1039.       K(98) = 1.
  1040.       K(99) = 10**(-11.23)
  1041.       K(100) = 10.**(-1.951)     
  1042.       K(101) = 5.0
  1043.       K(102) = 1.11D-11
  1044.       K(103) = 3.D-13
  1045.       K(104) = 1.2D-11
  1046.       K(105) = 0.85D-11
  1047.       K(106) = 7.7D-18
  1048.       K(107) = 1.7D2
  1049.       K(108) = 10.**(-10.62)
  1050.       K(109) = 10.**(-10.92)
  1051.       K(110) = 10.**(-12.22)
  1052.       K(111) = 2.6D-18
  1053.       K(112) = 3.3
  1054.       K(113) = 10.**(-11.23)
  1055.       K(114) = 10.**(-1.951)
  1056.       K(115) = 5.00135
  1057.       K(116) = 5.00135
  1058.       K(117) = 5.005
  1059.       K(118) = 10.**(-10.62)
  1060.       K(119) = 10.**(-12.1)
  1061.       K(120) = 10.**(-11.30)
  1062.       K(121) = 10.**(0.2553)
  1063.       K(122) = 3.0D-18
  1064.       K(123) = 3.0D-18
  1065.       RETURN
  1066.       END
  1067.  
  1068.  
  1069.       SUBROUTINE DEPOSNOP(VD,H,V)
  1070.       IMPLICIT DOUBLE PRECISION (A - Z)
  1071.       INTEGER I
  1072.       DOUBLE PRECISION VD(67),V(67)
  1073.       DOUBLE PRECISION J(79,57),ALFA
  1074.       COMMON /INTENS/ J,ALFA
  1075.  
  1076.       IF (ALFA.GE.86.) THEN
  1077.        V(1)=0.5*VD(1)
  1078.       ELSE
  1079.        V(1)=VD(1)
  1080.       ENDIF
  1081.       DO 314 I=2,67
  1082.       IF (ALFA.GE.86.) THEN
  1083.        V(I)=0.25*VD(I)
  1084.       ELSE
  1085.        V(I)=VD(I)
  1086.       ENDIF
  1087. 314   CONTINUE
  1088.  
  1089. C Deposition velocities divided by mixing height H => s-1
  1090.       DO 315 I=1,67
  1091.       V(I)=V(I)/H
  1092. 315   CONTINUE
  1093.       RETURN
  1094.       END
  1095.  
  1096.       SUBROUTINE EMISNOP(EMI1,H,EMI)
  1097.       IMPLICIT DOUBLE PRECISION (A - Z)
  1098.       INTEGER I
  1099.       DOUBLE PRECISION EMI1(67),EMI(67)
  1100.       DOUBLE PRECISION J(79,57),ALFA
  1101.       COMMON /INTENS/ J,ALFA
  1102.  
  1103.       DO 316 I=1,67
  1104.        IF (ALFA.LT.86.) THEN
  1105.     EMI(I)=1.5*EMI1(I)
  1106.        ELSE
  1107.     EMI(I)=0.5*EMI1(I)
  1108.     EMI(42)=0
  1109.        ENDIF
  1110. 316   CONTINUE
  1111.  
  1112. C Emission rates divided by mixing height H => molecule cm-3 s-1
  1113.       DO 317 I=1,67
  1114.     EMI(I)=EMI(I)/H
  1115. 317   CONTINUE
  1116.       RETURN
  1117.       END
  1118.  
  1119.       SUBROUTINE washout(T,H2O,WOUT)
  1120. c calculates the washout coefficient (s-1)
  1121.       implicit double precision (a - z)
  1122.       psat = DEXP(77.344913-7235.4247/T-8.2*DLOG(T)+0.0057113*T)
  1123.       rh=100.*H2O*T/psat*1.381d-17
  1124.       if (rh.gt.100.) rh=100.
  1125.       if (rh.ge.95.) then
  1126.     wout=1.d-4
  1127.       else 
  1128.     wout=0.
  1129.       endif
  1130. c      write(6,*) 'real RH= ',rh, 'wout= ',wout
  1131.       return
  1132.       END
  1133.  
  1134.       SUBROUTINE CLOUDINESS(CL,NHcl,ALFA,C1,koodi)
  1135.       IMPLICIT DOUBLE PRECISION (A-Z)
  1136.       integer cl,nhcl,koodi
  1137.       koodi=0
  1138.       if (alfa.gt.86.) then 
  1139.        alfa1=86.
  1140.       else
  1141.        alfa1=alfa
  1142.       endif
  1143.       M1=950./(1013.*COS(alfa1*3.14159265/180.))
  1144. c      write(6,*) 'm1= ',m1
  1145.       IF (CL.EQ.0) THEN
  1146.        C1=1.
  1147.       ELSE IF (CL.EQ.1) THEN
  1148.        T1=0.2684-0.0101*M1
  1149.        C1=1.-NHCL/8.*(1.-T1)
  1150.       ELSE IF (CL.EQ.2.OR.CL.EQ.3) THEN
  1151.        T1=0.3658-0.0149*M1
  1152.        C1=1.-NHCL/8.*(1.-T1)
  1153.       ELSE IF (CL.EQ.4) THEN
  1154.        T1=0.2363+0.0145*M1
  1155.        C1=1.-NHCL/8.*(1.-T1)
  1156.       ELSE IF (CL.EQ.5) THEN
  1157.        T1=0.4130-0.0014*M1
  1158.        C1=1.-NHCL/8.*(1.-T1)
  1159.       ELSE IF (CL.EQ.6) THEN
  1160.        T1=0.5456-0.0236*M1
  1161.        C1=1.-NHCL/8.*(1.-T1)
  1162.       ELSE IF (CL.EQ.7) THEN
  1163.        T1=0.8717-0.0179*M1
  1164.        C1=1.-NHCL/8.*(1.-T1)
  1165.       ELSE IF (CL.EQ.8) THEN
  1166.        T1=0.9055-0.0638*M1
  1167.        C1=1.-NHCL/8.*(1.-T1)
  1168.       ELSE 
  1169.        WRITE(6,*) 'TARKISTA SYÖTTÖARVOT!'
  1170.        koodi=1
  1171.       ENDIF
  1172.       RETURN
  1173.       END
  1174.  
  1175.       DOUBLE PRECISION FUNCTION TIMECUM(YEAR,iDAY,iMON,iHOUR)
  1176. c this function changes the date into the time calculated from the 
  1177. c beginning of the year in units "vuorokausi"     
  1178.       IMPLICIT DOUBLE PRECISION (A - Z)
  1179.       INTEGER YEAR,IDAY,IMON,IHOUR,ID,I
  1180.       INTEGER MON(12)
  1181.       DATA MON/0,31,28,31,30,31,30,31,31,30,31,30/
  1182.       
  1183.       ID=IDAY
  1184.       if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
  1185.      *mod(year,400).eq.0) MON(3)=29
  1186.       DO 172 I=1,IMON
  1187.       ID=ID+MON(I)
  1188. 172   CONTINUE
  1189.       TIMECUM=ID+IHOUR/24.                   
  1190.       RETURN
  1191.       END
  1192.  
  1193.       SUBROUTINE DATE(YEAR,TICU,PV,KU,HOUR)
  1194. C this function changes the time back to the form day month hour
  1195.       IMPLICIT DOUBLE PRECISION (A - Z)
  1196.       INTEGER I,YEAR,PV,KU,HOUR
  1197.       INTEGER MONTH(12)
  1198.       data MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
  1199.       
  1200.       I=1
  1201.       ID=0
  1202.       if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
  1203.      *mod(year,400).eq.0) MONTH(2)=29
  1204.  
  1205. 173   ID=ID+MONTH(I)
  1206.       IF (ID.LT.INT(TICU)) THEN
  1207.     I=I+1
  1208.     GOTO 173
  1209.       ENDIF
  1210.       PV=INT(TICU)-ID+MONTH(I)
  1211.       KU=I
  1212.       HOUR=(TICU-INT(TICU))*24+1
  1213.       RETURN 
  1214.       END
  1215.  
  1216.       SUBROUTINE SOLAR(RLAT,RLON,year,iDAY,iMON,HOUR)
  1217. C CALCULATION OF SOLAR FLUX S (W/M2), DECLINATION DEK (RAD) AND
  1218. C SOLAR HEIGHT ANGLE, GIVEN DATE,TIME (UTC),LAT(deg N),LONG (deg E):
  1219. C FORMULAE FROM PALTRIDGE AND PLATT (1977)
  1220.       IMPLICIT DOUBLE PRECISION (A-Z)
  1221.       INTEGER ALAT,ALON,year,IDAY,IMON,HOUR,ID,I,K,J
  1222.       INTEGER MON(12),month(12)
  1223.       DATA MON/0,31,28,31,30,31,30,31,31,30,31,30/
  1224.       ALAT=NINT(RLAT)
  1225.       ALON=NINT(RLON)
  1226.       data month/31,28,31,30,31,30,31,31,30,31,30,31/
  1227.       
  1228.       OPEN(11,FILE='ZENI.DAT', STATUS='UNKNOWN')
  1229.       FI=ALAT/57.3
  1230.       TICU=TIMECUM(YEAR,IDAY,IMON,HOUR)
  1231.       ticu=ticu-2./24.
  1232.       call date(YEAR,ticu,iday,imon,hour)
  1233.       ID=IDAY-1
  1234.       if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
  1235.      *mod(year,400).eq.0) MON(3)=29
  1236.       DO 170 I=1,IMON
  1237.       ID=ID+MON(I)
  1238. 170   CONTINUE
  1239.       if (mod(year,4).eq.0.and.mod(year,100).ne.0.or.
  1240.      *mod(year,400).eq.0) MONTH(2)=29
  1241.       DO 300 K=1,month(imon)
  1242.       D=ID*6.2832/365.
  1243.       DO 200 J=1,24-HOUR
  1244. C     S=1368.*(1.00011+.034221*COS(D)+.00128*SIN(D)+.000719*COS(2*D))
  1245.       DEK=.006918-.399912*COS(D)+.070257*SIN(D)-.006758*COS(2*D)
  1246.      C +.000907*SIN(2*D)-.002697*COS(3*D)+.00148*SIN(3*D)
  1247.       H=HOUR+ALON*24/360.
  1248.       A=COS(FI)*COS(DEK)*COS((H-12)*6.2832/24)+SIN(DEK)*SIN(FI)
  1249.       write(11,107) IDAY,IMON,HOUR,min(90.d0,90.d0-DASIN(A)*57.3)
  1250.  107  FORMAT(3I4,1F8.2)
  1251.       HOUR=HOUR+1
  1252. 200   CONTINUE
  1253.       iday=iday+1
  1254.       HOUR=0
  1255.       ID=ID+1
  1256. 300   CONTINUE
  1257.       CLOSE(11)
  1258.       END
  1259.